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ABSTRACT. In this paper, we discuss the application of quasi-Monte Carlo methods to 
the Heston model. We base our algorithms on the Broadie-Kaya algorithm, an exact sim- 
ulation scheme for the Heston model. As the joint transition densities are not available in 
closed-form, the Linear Transformation method due to Imai and Tan, a popular and widely 
applicable method to improve the effectiveness of quasi-Monte Carlo methods, cannot be 
employed in the context of path-dependent options when the underlying price process fol- 
lows the Heston model. Consequently, we tailor quasi-Monte Carlo methods directly to the 
Heston model. The contributions of the paper are threefold: We firstly show how to apply 
quasi-Monte Carlo methods in the context of the Heston model and the SVJ model, sec- 
ondly that quasi-Monte Carlo methods improve on Monte Carlo methods, and thirdly how 
to improve the effectiveness of quasi-Monte Carlo methods by using bridge constructions 
tailored to the Heston and SVJ models. Finally, we provide some extensions for computing 
greeks, barrier options, multidimensional and multi-asset pricing, and the 3/2 model. 



1. Introduction 

In this paper, we show how to apply quasi-Monte Carlo (QMC) methods to price 
path-dependent contingent claims where the underlying asset price is given by the Hes- 
ton model 1 28 1 and the SVJ model |8|. We recall that, on a filtered probability space 
(£2, J^, (j^" r ) r >o,P) under the assumption that P is (already) the risk-neutral pricing mea- 
sure, the Heston model is given by the system of stochastic differential equations 

dS, = rS,dt + VV,S,dW t l , 

r- 2 W 
dV t = k{9 - V,)dt + a Vv,dW r 2 , 

where (W t )t>o and (W, 2 ) r >o are Brownian motions under P with d(W l ,W 2 ) = pdt. The 
process S := {St)t>o models the asset price dynamics and V := (Vr)r>o the (stochastic) 
variance of S. Here, r is the risk-free rate of interest, is the long-term average variance, 
K is the mean-reversion speed of V, and <7 is the volatility of V. The SVJ model adds 
jumps in the dynamics of the underlying S and provides a model that calibrates better 
to the observed market prices of short-dated European options exhibiting steep skew, see 
Section [6] An important feature of our approach is that we can work in regimes where 
the Feller condition 2k9 > d 1 is violated, which is useful for market practitioners, see e.g. 

ED- 

The Heston model assumes diffusion dynamics for both the spot price and the volatility, 
but even though the spot price and volatility are jointly Markovian, we do not deal with 
a stationary, independent increment process. Therefore our situation is different to the 
majority of papers that apply QMC methods to finance problems 0] [2] [3] |4] [9] [14] |27j 
[3T1 [32l [38l |46l |49ll and, in particular, the Linear Transform (LT) method of Imai and Tan 
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OTl [32 1 (see also Leobacher ll38l ) does not seem applicable in our case of path-dependent 
options when the underlying follows either the Heston model or the SVJ model. 

We recall that the LT method relies on the observation that multiplying a vector of stan- 
dard normal random variates by an orthogonal matrix produces another vector of standard 
normal random variates. Consequently, evaluating the normal cumulative distribution func- 
tion elementwise at the vector resulting from the multiplication of an orthogonal matrix by 
a vector of standard normal random variates, produces a vector of random variates uni- 
formly distributed on [0, 1]. One consequently uses these uniform random variates in the 
QMC procedure relevant to the problem. The key step of the LT method is a judicious 
choice of the orthogonal matrix. In fact, an optimization procedure is used to obtain this 
matrix, see OTl . 11321 . Unfortunately, in [32] the method is only presented for independent 
increment processes (see Proposition 4. 1 in 1321 ) and requires knowledge of the transition 
density of the underlying stochastic processes. Though the density of the spot price, con- 
ditional on the initial values of spot and volatility, in the Heston model has been studied 
1 7 1, the joint transition density of spot and volatility, which is needed for the valuation of 
path-dependent options, does not seem to be available. Consequently, the LT method is not 
applicable in our situation, instead we directly construct bridges for the stock price process 
under consideration. 

We recall that QMC is a class of numerical methods for high-dimensional integrals that 
can broadly be divided into two categories: nets [21,42] and lattice rules [42,48]. However, 
in practice, once the problem is properly formulated both approaches may be applied. This 
is the first of three contributions of the paper: we first show how to formulate the finance 
problem as an integration problem based on the exact simulation scheme of Broadie and 
Kay a [ 12] then, similar to the results from the references cited above, we demonstrate that 
QMC methods outperform Monte Carlo (MC) methods when using Sobol point sets with 
Owen's scrambling method. We also improve our QMC results slightly by conditioning in 
the case of European call options. 

Secondly, we extend the results for the Heston and SVJ model to the case of path- 
dependent options and demonstrate that by allocating more of the variance to the early 
dimensions of the QMC point set, the performance of QMC methods can be significantly 
improved. This is also in line with the results presented in the cited references. We achieve 
this by employing a bridge construction for the variance process based on ll40l I5T1 and 
the well-known bridge construction for a Brownian motion with drift and time-dependent 
volatility 11241 and the bridge construction for jump processes [5]. These constructions 
are in the spirit of l2l [51 l46l . where bridge constructions were derived for the stochastic 
processes under consideration. We demonstrate our bridge algorithm by considering Asian 
call options in the case of the Heston and SVJ model. 

Thirdly, we provide a number of additional results of interest in finance. We show 
how to compute greeks in our framework, we consider Barrier options using the ideas of 
Glasserman and Staum ll26l . multi-asset stochastic volatility models, and the extension to 
the 3/2 model. 

The remainder of the paper is structured as follows. In Section|2j we discuss the Broadie- 
Kaya algorithm, which forms the basis of the algorithms introduced in this paper, and 
show how to combine it with the QMC methods. In Section [3] we recall QMC methods 
and Section |4] shows how to effectively combine the Broadie-Kaya algorithm with QMC 
methods in the context of path-independent European options. Section[5]introduces bridge 
sampling for the square-root process and the spot price process. The extension to the SVJ 
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model is discussed in Section [6] Further extensions of the algorithms presented in this 
paper are discussed in Section|7] 

2. QUANTILE FUNCTION FOR THE HESTON MODEL 

To price any contingent claim using a MC or QMC method, one must be able to sample 
from the law of the asset price process in a fast and accurate way. Our approach is to ap- 
ply the inverse transform method which requires us to obtain the quantile function of the 
law of the Heston process (S t )t>o- We recall that given a random variable X with distribu- 
tion function F(x) := P(X < x), the quantile function (or inverse cumulative distribution 
function) returns the value below which random samples from the given distribution would 
occur p amount of the time: the quantile function Q : [0, 1] — > K for the distribution of X is 
defined as Qx(p) := inf{x eR:p< F(x)}. 

For t > u, the distribution of St and V, solving (Q]) with initial conditions S u E R+ and 
V u € K + can easily be shown to be given by 

S t = S u exp (r{t - u) - - jf ' V s ds + p jf W<dW; + ^/l-p 2 £ \/V~ s dW s l 

V, = V u + K6{t- s) - kJ V u du + a J V%dW^. (2) 

Our approach to obtain the quantile function for the distribution of the random variable S t 
is based on the exact simulation method obtained by Broadie and Kaya fl2l . 

2.1. The Broadie-Kaya approach. We recall that the exact simulation approach for (Q]) 
given by Broadie and Kaya lfl2ll is as follows: 

(1) Simulate V t given V u = x, 

(2) Generate a sample from the distribution of J" V s ds given V u = x and V t = y, 

(3) Recover j' u ^/V~ s dW} from (f2]l given V, = y, V u = x, and J u V s ds — zas 

f Vv7,dW* = -(y-x+Kd(t-u)-z), 

(4) Generate a sample from the distribution of S, given f u ^/V s dWg and J' u V s ds. 
When simulating a random variable X using QMC methods it is sometimes more con- 
venient to rephrase such algorithm in terms of quantiles, allowing one to substitute either 
random uniforms (for a MC method) or QMC point sets in the quantile function Qx to 
draw a sample of X. 

2.2. Obtaining quantile functions. The Broadie-Kaya approach can be reformulated in 
terms of quantiles. To simplify notation, we shall henceforth use the following convention: 
we write X ~ L to denote that X has a law L and X\Y —y^Lto denote that the random 
variable X conditioned on the event Y — y has law L. We sometimes write X \y instead of 
X\Y = y when the context is clear. 

The quantile function for V t \V u — x, which is needed for step 1 in the Broadie Kaya 
approach, is easily found. It is well known that for t > u the distribution of V t \V u follows a 
noncentral chi-squared distribution: 

VIV - r- g2 ( 1 - eX P("' C ( f - M ))) r 2 

v t \v u -x — Xd 

where Xv(^) denotes the noncentral chi-squared random variable with v degrees of free- 
dom and noncentrality parameter A. Therefore, the quantile function for V t \V u = x can 



4fcexp(— K(t — u)) 
<7 2 (1 -exp(— jf(f — «)))" 



d:= 



49k 
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be obtained from the quantile function for a noncentral chi-squared distribution with the 
appropriate choice of parameters V and A . 

The quantile function of f V s ds\Vu = x, Vt = y is the most expensive step. Broadie and 
Kaya obtained its distribution function by computing the conditional Laplace transform 
and performing a numerical inversion of the characteristic function. An alternative ap- 
proach to determining the distribution function was given by Glasserman and Kim using a 
series expansion of the random variable [25]. We follow the approach given by Broadie and 
Kaya whereby the conditional distribution F(x) := P ( / M V s ds < x\V u ,V t ) is obtained by nu- 
merically inverting the (conditional) characteristic function <$>(a) :— E [exp (ia j u V s ds \V u ,Vt)] . 
We recall that <t> was given explicitly as 

4> (a) = — ri-, — rr --^exp C(a 

v ' k(i _ e -rW('-»)) B K " 

where the terms A, B, C, and 7(a) are 




(7 2 (1 - e -rW('-«)) 

AKe \K(t-u) \ 



C{a) : 



. e -<<-")) / ' 

V u + V, ( _ 7(a)(l+ e -yM('-")) ' 

a 2 \ l-e- K (*-«) l-er(f)(<-«) 



7(0) := \/ k 2 - 2a 2 ia, 

and I v (x) denotes the modified Bessel function of the first kind. The inversion of <f> (which 
gives F(x)), can be calculated by applying a trapezoidal rule to approximate the integral 

F(*)=!r^H»fr)A, 

it Jo z 

where 9\(x + iy) := x for x,y £ R. We note that the term B does not depend on a so only 
needs to be computed once per inversion of <t>. Given a mesh size of h, the integral is 
approximated as 

F(x) = hjc + 2 ^ 
% %p x j 

where £/, is the discretization error associated with the choice of mesh size h and £n is 
the truncation error caused by taking N <°°. Following Broadie and Kaya, to approximate 

F(x) with good accuracy we seth = 2n/ (x+\m\ \+q\^Jni2 — m 2 \) where m.\ and 1112 are the 

first and second moments of J u V s ds\V u ,Vt and q £ N is chosen sufficiently large (e.g., q = 
5). Explicit (but long) expressions for m\ and mi can be found using a computer algebra 
system (CAS) using the well-known technique of differentiating <t> and setting a = 0. The 
upper bound of the summation is chosen to satisfy \<P(hN)\/N < ne/2 where e is the 
desired truncation error. Thus, given our numerical approximation F h>N of the distribution 
function F we can apply a root finding procedure to identify the quantile function of F h ' N . 

Finally, we consider the quantile function of S t conditioned on the event f u V s ds = y, the 
event J u ^/V^dW 2 = z, and the initial condition S u . First, observe that yl — p 2 ./,' VVsdW s l 
is normally distributed with mean zero and variance f u V s ds. Therefore, \og(S t )\y 7 z 7 S„ 
is normally distributed with mean r(t — u) — jy + pz and variance (1 — p 2 ) f u V s ds — 
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(1 -p 2 )y- Further, if we have the quantile function gz forZ~ N(0, 1) then if X ~ N(fl, a 2 ) 
then Qx{p) = M + <*Qz(p) and Q exp(X ) (p) = exp(ju + <7Qz(p)) as e* is an increasing func- 
tion of x. Hence, we have the quantile function for S t \y,z,S u . 



3. QMC METHODS 

Regarding QMC point sets, we will make use of digital nets II2T1 l42l . Given a (base 
b) digital net (u,)" =1 , where u, £ [0, l] d for some dimension d £ Z + , we will always ran- 
domize the point set using Owen's scrambling algorithm to compute standard errors P3l . 
Given a generic point x £ [0, \ ) d , where x = (x\,... ,x s ), we recall that Owen's algorithm 
expands each Xj as 

ft + ft 2 + ---' 

and generates a scrambled point y £ [0,1 ) , where y = (y\ , . . . , yj), as 

»-^ + - 

The permutation 7T, applied to j, j = l,...,d depends on for 1 < /. In particular, 
'?/.: = w>(§y,i), »?,.; = T^fe), iJ/,3 = ^^l^C^'-s) and in general 

where 7tj and 7T ; - ^ , ^ t _ ( , ^ > 2, are random permutations of {0, 1 , . . . , b — 1 }. We assume 
that permutations with different indices are mutually independent. We recall that if the 
scrambling algorithm is applied to x to obtain y, then y is uniformly distributed in [0, l) d 
by Proposition 2 in |43 1. We find it convenient to introduce the following notation 

y = ^(x) , 

so n(-) represents the scrambling algorithm applied to x to obtain y. 

We recall that by our reformulation of the Broadie Kaya approach for the Heston model 
in Section [2] in terms of quantiles, the discounted payoff e~ rT g(Sj) for some T > can 
be rewritten as a function / : [0, l] 3 — > K and, in the case of a path-dependent payoff, 
e~ rT g(S tl ,S t2 ,..., S t „ ) for times < fi , t2, ■ ■ ■ , t„ = T can be transformed to the function 
/ : [0, l] 3x " — > R. Henceforth, we simply assume that all discounted payoffs are mapped 
to the cZ-dimensional (for some appropriate d) unit cube in this manner. 

Given q independent permutations n r , r = l,...,q, and the discounted payoff of the 
financial derivative / : [0, \] d — > R, we estimate the price of the derivative using 

1 1 1 1 1 n 

/QMC = -E/r=-E-E/(^(u < )), 

and compute standard errors via 



0"QMC 



q(q-l) 



We point out that for digital nets, scrambling is the preferred randomization method, as it 
achieves the optimal convergence rate J6] |2T] |45) . For an implementation of the scrambling 
algorithm, we refer the reader to ||30l . 

For purposes of comparison, we will also look at MC estimators. In this case, we will 
choose q x n points, (t^i)'-*", independent and identically distributed in [0, l] d , and estimate 
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derivative prices using 



-. qxn 

toe = — !/(& 



and compute standard errors using 



0MC — 



'Ef=T(/te-)-toc) 2 



q x «(g xn-1) 

We conclude the section by commenting on the variances of /mc an d ^qmc- For square- 
integrable functions /, it is well-known that 

Var(/ M c) = • 

qxn 

Regarding Iqmc, it is known that for a square-integrable function / the Monte Carlo conver- 
gence rate is matched P41I2T1 . However, scrambled digital nets can exploit the smoothness 
of the integrand: if the integrand satisfies a Holder condition of order a with < a < 1, 
then the variance decays at a rate of order ff(n~( l+2a " ,+E ) for some e > 0, and the leading 
constant is allowed to depend on the dimension of the point set [6 .21]. In practise, it is dif- 
ficult to confirm the smoothness of the integrand under consideration, hence it is important 
to investigate the standard errors of QMC methods numerically. 

4. European Call options 

In this section, we follow [50] and apply conditional Monte Carlo method to improve 
the efficiency of the Broadie-Kaya algorithm. We use European call options to demonstrate 
the method. The approach relies on the observation that given 

rT rT 

/ V s ds, / WsdW? 
Jo Jo 

the price of a European call option is given by the Black-Scholes price, with modified 
initial share price 

So = S exp(-Y J T V s ds + pJ T y/VsdW* 



and adjusted volatility d\J 1 — p 2 , where a := y j J V s ds. Using BS(So,K,r,T,a) to 
denote the Black-Scholes price of a European call with initial stock price So, strike K, 
interest rate r, time to maturity T and volatility a , we have 

E(e- rT {S T -K) + ) =E(BS(So,K,r,T,6)) , 

that is, we firstly simulate J T V s ds and J T \/VsdW} using the Broadie-Kaya algorithm, and 
then compute the Black-Scholes price, for the particular values of So and a corresponding 
to J T V s ds and J r \/V s dW}. In Table [T] we show price and standard error estimates for 
a European call option. Columns 1 and 2 show that combining QMC methods with the 
Broadie-Kaya algorithm already improves on the Monte Carlo method. However, using 
the conditioning argument, the QMC method significantly improves on a naive application 
of QMC methods to the Broadie-Kaya algorithm. There are two reasons for the variance re- 
duction: we do not estimate the Black-Scholes price using Monte Carlo simulation, which 
is done in Step 4 of the Broadie-Kaya algorithm, but compute the value exactly. Further- 
more, it can be expected that when combining the conditional expectations with QMC 
methods, the approach is even more efficient [50]. This is due to the fact that taking the 



QUASI-MONTE CARLO METHODS FOR THE HESTON MODEL 



7 



Table 1. Estimated prices (with standard errors given in parentheses) 
of European options where the asset price process is given by a Heston 
model. These values are based on 30 independent batches. 



Trials 



MC 



QMC 



Cond QMC 



128 6.847642(0.121792) 

256 6.785936 (0.085627) 

512 6.794658(0.059212) 

1024 6.818196(0.042269) 

2048 6.820823 (0.030005) 

4096 6.815857(0.021199) 

8192 6.789039(0.014945) 

16384 6.800179(0.010576) 



6.792353 (0.074884) 
6.815575 (0.037468) 
6.807123 (0.023326) 
6.805009 (0.010747) 
6.805464 (0.004077) 
6.806346 (0.002346) 
6.806315(0.001721) 
6.806484 (0.000730) 



6.807731 (0.011452) 
6.807578 (0.007037) 
6.806918(0.001924) 
6.807080(0.001928) 
6.806219(0.001182) 
6.806326 (0.000480) 
6.806558 (0.000520) 
6.806438 (0.000249) 



Option parameters: 5 = 100, K = 100, V 
r = 3.19%, T= 1.0 



= 0.010201, K = 6.21, e 
year, true option price is 



= 0.019, a = 0.61, p = -0.70, 
6.80611. 



conditional expectation has a smoothing effect, which can be expected to improve the per- 
formance of quasi-Monte Carlo methods, see 11371 Subsection 10.1]. 

5. Bridge sampling for path-dependent options 
In this section, we study the pricing of options whose payoff is a function of 

Sfl 5 • • • ) St h , 

where we choose h = 2"' for simplicity. A naive approach to simulating this path is given 
in AlgorithmQ] The algorithms for the simulation of the random variables required in steps 
1, 2, and 3 of Algorithm Q] (using either MC or QMC) follows from the algorithm given in 
Section|2] 

Algorithm 1 Naive version of the Broadie-Kaya exact simulation algorithm 
1: Simulate (V ti )f =1 in the following order 

V tl ,V h ,---,V tN . 

2: Simulate (j^ V s ds\V, i _ l , V f ,.) for i = 1,...,N. 
3: Simulate (S ti )f = i in the following order 

St 1 , St 2 , St 3 , • • • , St N . 



It is well-known that QMC methods are particularly efficient if the effective dimen- 
sion of the problem under consideration is low 1141 . In our case, the effective dimension 
depends on the variance of the discounted pay-off of the financial derivative under consid- 
eration. To reduce the effective dimension, following ll2ll5l fT4l . we now allocate the early 
dimensions to variates with high variances. We achieve this by proposing a bridge algo- 
rithm, given by Algorithm^ to generate the paths of (S t ) t >o at the time points t\ , . . . ,t^- 
The simulation of the random variables required in our bridge algorithm requires the con- 
struction of a bridge sampling algorithm for the square-root process (V^),>o and a bridge 
sampling algorithm for the stock price process (S t ) t >o- In the coming sections, we propose 
algorithms to construct these bridges. 
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Algorithm 2 Bridge version of the exact simulation algorithm 



1: Simulate (Vi,)£i in the following order 

V, V, V, V, 



2: Simulate (j^ V s ds\V ti _ n V t ,j for z = 1, . . . ,N. 
3: Simulate {S tj )f =1 in the following order 



5.1. Bridge sampling for square-root processes. We now recall how to perform bridge 
sampling for square-root processes by relying on the bridge sampling algorithm for squared 
Bessel processes studied by Yuan and Kalbfieisch [51 1 and Makarov and Glew ll40l . The 
following well-known result linking square-root processes and squared Bessel processes is 
employed, see Proposition 6.3.1.1 in ll34l . 

Proposition 1. The square-root process V — {V t , t > 0} is a squared Bessel process trans- 
formed by the following space-time change: 



where c(t) = j^(e KI — 1), X = {X t , t > 0} is a squared Bessel process of dimension 8 = 
4k9 



As we assume that K, 9, a > 0, we have that for 8 > 2, the boundary is polar and for 
< 8 < 2 it is reflecting, see Figure 6.1 in ll34ll . This allows us to propose Algorithm [3] 
adapted from the algorithm in [40], to sample a square-root process at times t\ , . . . , f/, where 
h = 2 m withm> 1. 

5.2. Bridge sampling for the stock price. We now propose a bridge sampling algorithm 
for step 3 of Algorithm [2] that is, an algorithm to simulate (5(, )| ! =1 where h = 2 m in the 
order S th , St h , 2 , St h , 4 , St 3h i 4 with effect of allocating the early dimensions to variates with 
high variances. The following lemma allows us to set up a sampling scheme for the stock 
price process. 

Lemma 1. Given times u <t\ <t^< ?3, the joint distribution o/log(5, 1 ),log(5'( 2 ),log(5', 3 ) 
conditional on the initial condition S u , the variances V :— (V tl , V t2 , V f3 ), and the integrated 
variances IV := (JJ V s ds, f'^V s ds, J u 3 V s ds) is given by 



V t =e- Kt X c{t) , 



.2 




wherefor i — 1,2,3 we have defined of := (1 — p 2 ) ffj V s ds and 



1 /•'< 

IM := log(S f0 )+rti-- I V s ds + p 



£ Jit Jit 




V s dW s \ 



Further, we have \og(S tl )\S tl ,S t3 ,V,TV ~ N(M,'L) where 




log(S n ) + Hi - Hi) 
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Algorithm 3 Bridge sampling for the square-root process 

Require: Time indices (t\ , . . . ,?2 m ) 
1: h <- 2"\ ; max <- 1, x <- V , 8 <- ^ 

2: for z = 1,.. . ,/i do sj «- f^(exp(K - f ! ) — 1) 

3: Generate Z ^ ^ 2 (5,xo/iAf) 

4: X/, ■<— SjjZ 
5: f <-0 

6: for k — 1 , . . . , m do 

7: i'min ^ / ^ z'min 

8: / «- 0, r <r- h 

9: for; = l,...,; max do 

1 1 : Generate P ~ Poisson (A ) 

12: Generate Z ~ Besselff - 1, 

13: Generate G ~ Gamma(P + 2Z + f , 2(^-^7(7,—^-) ) 

14: X; «- G 

15: i <- i + h, I <— l+h, r <— r + h 

16: ./max ^ 2y'm ax , /l < z'min 

17: for / = 1 , . . . , h do Vi '<— exp(— KtijXi 

return sampled path (vi , . . . , v/,) from distribution of (V f] , . . . , V tim ) 



(a 2 (t 2 )-a 2 (tt)f 



E:=a z (f 2 )-ff 2 (fi)- 



a 2 (f 3 )~cr 2 (fi) 



Proof. The first part of the lemma is straightforward, the second follows immediately from 
the conditioning formula (2.25) in [24 1. □ 

By applying Lemma[T] we propose Algorithm|4]to sample the path of the share price pro- 
cess (S t ) t >Q at the time points (S tl ,...,S th ). This algorithm requires the drifts (fli , . . . , )J.2 m ) 
and volatilities (of, . . . , ofm) determined in LemmaQ] 

5.3. Numerical results. To demonstrate our bridging technique, we estimate the price of 
an Asian call option where the asset price process is given by a Heston model. We recall 
that the payoff of an Asian call option, for strike price K and expiry T, is given by 




where ?i,...,f^ are d monitoring dates. We compare the prices obtained with standard 
Monte Carlo, QMC using Sobol points and Owen's scrambling algorithm, and the bridge 
construction using Sobol points and Owen's scrambling. These results are presented in 
Table |2] We find that QMC methods already improve on Monte Carlo methods, however, 
using bridge constructions, the effectiveness of QMC methods can be enhanced. 

6. Extension to the SVJ model 



It is well known that Heston stochastic volatility model cannot fit short-term smiles well 
if they exhibit skew [23 1. This motivates the introduction of jumps into the dynamics of 
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Algorithm 4 Bridge sampling for the share price process 

Require: Time indices (t\,.. . ,t%™), drifts (jUi,-- .,jJ.2"<), volatilities (of, . . .,0^), and So 
l:h<- 2 m , ; max <- 1 
2: Generate (Z h . . . ,Z h ) ~ 2V(0,7) 

3: S A «- exp(jU, 7 + OfcZ A ) 
4: S <- ^0 

5: for k= 1, . . . ,m do 

6: i min <- ft/2, / «- i min 

7: Z <- 0, r <- ft 

8: fory = l,...,; max do 

2_ 2 

9: fl <- jll; - JU/ + logfo) + 4— 4 (log(s r ) - log(if/) Mr) 

°r °; 

10: b<-Of-Of- S 2 

' ' Gf— of 

11: Si <— exp(a + bZj) 

12: / <— / + ft, / -s— / + ft, r -s— r + ft 

13: 7max ^ 2y max , ft < 'min 

return sampled path (si,..., S21" ) from distribution of (5 fl , . . . , 5( 2m ) 



Table 2. Standard errors of Asian option price where the asset price 
process is given by a Heston model using Monte Carlo (M), Quasi-Monte 
Carlo (Q), and the Bridge Quasi-Monte Carlo (B) algorithms. These 
values are based on 30 independent batches and m trials. 



n 


m: 


2 b 


2> 


2 B 


2 9 


2 1U 


2 11 


2 12 


2 1J 


2 14 




M 


0.0959 


0.0740 


0.0639 


0.0365 


0.0302 


0.0204 


0.0128 


0.0086 


0.0052 


4 


Q 


0.0303 


0.0231 


0.0147 


0.0079 


0.0039 


0.0033 


0.0020 


0.0009 


0.0008 




B 


0.0299 


0.0172 


0.0133 


0.0070 


0.0053 


0.0036 


0.0028 


0.0025 


0.0010 




M 


0.1054 


0.0713 


0.0399 


0.0338 


0.0253 


0.0166 


0.0126 


0.0075 


0.0052 


8 


Q 


0.0413 


0.0285 


0.0160 


0.0080 


0.0056 


0.0043 


0.0023 


0.0021 


0.0009 




B 


0.0379 


0.0280 


0.0149 


0.0151 


0.0063 


0.0032 


0.0026 


0.0021 


0.0011 




M 


0.0973 


0.0646 


0.0435 


0.0351 


0.0244 


0.0227 


0.0127 


0.0074 


0.0069 


16 


Q 


0.0462 


0.0241 


0.0189 


0.0102 


0.0073 


0.0052 


0.0027 


0.0019 


0.0010 




B 


0.0409 


0.0298 


0.0182 


0.0130 


0.0061 


0.0042 


0.0029 


0.0018 


0.0013 




M 


0.0804 


0.0608 


0.0432 


0.0364 


0.0269 


0.0163 


0.0106 


0.0078 


0.0055 


32 


Q 


0.0523 


0.0377 


0.0181 


0.0185 


0.0123 


0.0095 


0.0035 


0.0022 


0.0018 




B 


0.0339 


0.0269 


0.0283 


0.0114 


0.0101 


0.0047 


0.0030 


0.0026 


0.0009 




M 


0.0913 


0.0656 


0.0362 


0.0349 


0.0246 


0.0158 


0.0127 


0.0080 


0.0064 


64 


Q 


0.0575 


0.0400 


0.0251 


0.0219 


0.0182 


0.0113 


0.0064 


0.0034 


0.0020 




B 


0.0537 


0.0205 


0.0118 


0.0069 


0.0076 


0.0061 


0.0046 


0.0033 


0.0039 




M 


0.0920 


0.0612 


0.0318 


0.0314 


0.0253 


0.0142 


0.0113 


0.0099 


0.0064 


128 


Q 


0.0705 


0.0331 


0.0319 


0.0208 


0.0201 


0.0123 


0.0054 


0.0038 


0.0027 




B 


0.0641 


0.0672 


0.0587 


0.0110 


0.0075 


0.0066 


0.0047 


0.0037 


0.0034 


Option parameters: S 


= 100, K 


= 100, Vb 


= 0.010201, K = 


6.21, e = 


0.019, <t 


= 0.61,p 


= -0.70, 



r = 3.19%, with n time monitors over the time period [0, 1]. 



the underlying share price process. The following model, often called the SVJ model, was 
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first proposed in 1 8 1 : 

dS t = S,- ((r-Xp,)dt + V t (pdW? + yjl-p z dwf\ + (Y, - l)dN^j , (3) 

dV, = k(Q- V,)dt + aVVtdW, 1 , 

where N, is a Poisson process with constant intensity X, (W t l ) r >o and (W r 2 ) r >o are indepen- 
dent Brownian motions both independent from the Poisson process (N t ) t >o, and the jump 
variables Y t are a family of independent random variables all having the same lognormal 
distribution with mean jis and variance a|. Furthermore, E(l* — 1) = p, so it follows that 

1 , 

Ai s = log(l+£)--cr s z . 

Solving the SDE for the stock price process (O, we obtain 

Nt 

St=S t Y{Yj, (4) 

where 

S t = 5 exp Ur - Xp)t - l - jf V s ds + p £ VV s dW} + W,dW? 

6.1. Exact simulation for the SVJ model. As discussed in 02|[35], equation (01 moti- 
vates the simulation algorithm for the SVJ model: First simulate the diffusion part as in 
Section|2]and then take care of the jump part using YljLi This results in Algorithm[5] 
which is the analogue of the Broadie Kaya algorithm presented in Section [2] for the SVJ 
model. We recall that this algorithm also appeared in |[T2l and in similar form as Algorithm 
7.1 in [ 35 1 . Further, since follow the lognormal distribution with mean n$ and variance 
a|, it is clear that 

Nt 

£log(tj)\N t ~N(NtH„N t o$) . 



Algorithm 5 Exact simulation algorithm for the SVJ model 

1: Simulate V t given Vq 

2: Generate a sample from the distribution of /j, V s ds given V t and Vq 

3: Recover y/V s dW} from © given V t , V and J ' V s ds 

4: Generate S t 

5: Generate Nt 

6: Generate Tjjii Yj, given N t 



There are alternative approaches to simulating IT/li Yj as required in Algorithm |5] As 
in Section 3.5 of J24), one can simulate N t by simulating the jump times of the Poisson 
process. Furthermore, as discussed in [ 12 1, given A',, one can simulate the jump sizes %, i = 
\....,N t individually. However, Algorithm|5]results in a problem that is of fixed dimension, 
in particular, the dimension of the problem in Algorithm[5]is five, i.e. five random uniforms 
(or a QMC point from [0, l] 5 ) are used to obtain a realization of St- Having a problem of 
fixed dimensionality is important when applying QMC methods, which is the ultimate goal 
of this paper, hence we choose the formulation presented in Algorithm!?] 
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6.2. Path-dependent options in the SVJ model. As in Section [5] we now turn to the 
problem of studying the pricing of options whose payoff is a function of 

where we choose h = 2™ for simplicity. A naive approach, analogous to Algorithm Q] is 
given by Algorithm^ 

Algorithm 6 Naive version of the exact simulation algorithm for the SVJ model 
1: Simulate (Vr,)f = i in the following order 

v h ,v h ,---,v th . 

2: Simulate ( V !i ds\V, i _ l , VU ,i=l,...,h. 

N 
i=l 



3: Simulate (^)^Lj in the following order 



4: Simulate (M, )f = i in the following order 



N tl ,N t2 ,...,N th . 



5: Simulate Ml^L^ +1 TylM,^ ,i = 1, ••• ,h, in the following order 



N n N, 2 



m> n ^ n 



6.3. Bridge sampling for path-dependent options in the SVJ model. Now, similar to 
the case of the Heston model presented in Section [5] we now propose an algorithm that 
allocates the early dimensions to variates with high variances. As such, we propose a 
bridge sampling algorithm for the SVJ model. We use the approach from Section|5]to deal 
with the diffusion component (§ tj )f =l and the approach proposed by Baldeaux in 13 to 
deal with the jump part. This results in Algorithm[7] 

The next lemma, also presented as Lemma 3.1 in @, allows us to perform Step 4 of 
Algorithm|7]efficiently. 

Lemma 2. Let s < u < t and k\ < k^. Then conditional on N s = k\ and N t = k2 the 
increment N„ — N s has the binomial distribution with parameters k^~k\ and (u — s)/(t — s). 

Finally, the next lemma provides the tool required to complete the final step of Algo- 
rithm [7] 

Lemma 3. Given times t\ < t% < ?3, the joint distribution of 

( N 'l N >2 N >3 \ 

£tog(?y),£log(Fy),£log&) 



conditional on (N tl ,N t2 ,N t3 ), is given by 



N,. ,N t2 ,Nh 



■ N 



Mi 

J"3 
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where /i, = N tj Hs an d of = Nt f Og. Further, we have that 



where 



and 



I>g(?;) 

7=1 



N n N, 3 

£l0g(fy),£l0g(fy)~^(M,E) 

7=1 ;=i 



°2 _ / x 

M := x\ +1112 — nil H — % k \X"i — x\ +m\ — m^) , 

o'-o{ 



2 _2 (g|- g l 2 ) 2 



£ := O2 — of — 



°3 - °1 



Algorithm 7 Bridge version of the exact simulation algorithm for the SVJ model 
1: Simulate (V^)^.j in the following order 

V, V, V, V, 

2: Simulate V^V^, V f j) , / = 1, . . . ,N. 
3: Simulate {S tj )f =l in the following order 

4: Simulate in the following order 

Nto,N tN ,N tN ,N t3N ,... 

5: Simulate (Ih= 1 ^/)jLi conditional on (N tj )f =l in the following order 



;=i ;=i .7=1 ;=i 



6.4. Numerical results. To demonstrate our bridging technique, we estimate the price of 
an Asian call option where the asset price process is given by a SVJ model. We com- 
pare the prices obtained with standard Monte Carlo, QMC using Sobol points and Owen's 
scrambling algorithm, and the bridge construction using Sobol points and Owen's scram- 
bling. These results are presented in Tableland agree with the ones presented in Table 
12 QMC methods improve on Monte Carlo methods, but the effectiveness of Monte Carlo 
methods can again be enhanced via bridge constructions. 

7. Further Extensions 

We now propose a number of extensions to the results of the previous sections, in par- 
ticular, the algorithms discussed in Section and Section |6]can be modified to solve some 
further problems of interest to practitioners in finance. First, we discuss how to compute 
"greeks". Second, we show how our algorithm can be enhanced for barrier options. Third, 
we provide an extension to the multidimensional and multi-asset setting. Fourth, we show 
how to price path-dependent options when the asset price process follows the 3/2 model. 
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Table 3. Standard errors of Asian option price where the asset price 
process is given by a SVJ model using Monte Carlo (M), Quasi-Monte 
Carlo (Q), and the Bridge Quasi-Monte Carlo (B) algorithms. These 
values are based on 30 independent batches and m trials. 



n 


m: 


2 6 


2> 


2* 


2 y 


2 lu 


2 n 


2 [1 


2 [i 


2 14 




M 


0.1354 


0.0904 


0.0717 


0.0525 


0.0254 


0.0199 


0.0152 


0.0105 


0.0065 


4 


Q 


0.0549 


0.0401 


0.0293 


0.0116 


0.0064 


0.0055 


0.0042 


0.0026 


0.0018 




B 


0.0595 


0.0403 


0.0217 


0.0214 


0.0171 


0.0041 


0.0032 


0.0029 


0.0015 




M 


0.1313 


0.0873 


0.0442 


0.0421 


0.0284 


0.0200 


0.0089 


0.0078 


0.0060 


8 


Q 


0.0488 


0.0308 


0.0304 


0.0198 


0.0188 


0.0089 


0.0031 


0.0027 


0.0016 




B 


0.0455 


0.0235 


0.0224 


0.0151 


0.0088 


0.0087 


0.0044 


0.0036 


0.0027 




M 


0.0786 


0.0708 


0.0447 


0.0427 


0.0278 


0.0199 


0.0127 


0.0107 


0.0060 


16 


Q 


0.0618 


0.0415 


0.0312 


0.0185 


0.0109 


0.0088 


0.0062 


0.0035 


0.0021 




B 


0.0520 


0.0349 


0.0196 


0.0173 


0.0100 


0.0069 


0.0058 


0.0023 


0.0015 




M 


0.1218 


0.0622 


0.0508 


0.0346 


0.0265 


0.0172 


0.0116 


0.0082 


0.0060 


32 


Q 


0.0615 


0.0318 


0.0353 


0.0213 


0.0146 


0.0084 


0.0046 


0.0037 


0.0020 




B 


0.0359 


0.0327 


0.0193 


0.0140 


0.0105 


0.0052 


0.0031 


0.0019 


0.0012 



Option parameters: 5 = 100, K = 100, V = 0.010201, K = 6.21, 6 = 0.019, a = 0.61, p = -0.70, 
r = 3.19%, X = 0.11, fa = —0.1391, and a s = 0.15, with n time monitors over the time period [0, 1]. 



7.1. Computation of greeks. In a subsequent paper by Broadie and Kaya [O), it was 
shown how to modify their exact simulation algorithm to compute greeks. We now sum- 
marize their methodology and show how to adapt their method to our bridge algorithms 
for path-dependent options. To simplify the exposition, we focus on the Heston case, but 
we can also handle the SVJ case. We discuss the pathwise (PW) and the Likelihood Ratio 
(LR) method ifTTl l24l . Using the notation from |T3j, we assume that the option price is 
given by 

a(d)=E\f(6)], 

where / denotes the discounted payoff function and 9 the parameter of interest, i.e., we 
are interested in computing a'(0). Regarding the PW method, we have 

a'(e) = -^E\f(e)}=E[f(9)] , 

assuming the interchange of differentiation and integration is permitted. For the LR method, 
we consider 9 as a parameter of the transition density of the random variable under consid- 
eration, say X. Denoting this density by ge(x), we have 

«'(a) = A £L/(z)] = ^ /(x) ^ e(x) , x . 

Now, we rewrite this as 

a'(9)= [ d m 8 -^ 
Jm. d ge (x) 

The quantity is also known as the score function, and is of course independent of the 
particular payoff under consideration. Clearly, both approaches rely on the interchange- 
ability of differentiation and expectation, and we refer the reader to IfTTl |24| for details. 
To be able to apply the LR method, we need to have access to the density, denoted by 
ge(x) above. To achieve this, we apply a conditioning argument, from the law of iterated 
expectations, 
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E\f(X)]=E[E[f(X)\Y}}, (5) 

where Y is a vector valued random variable. For the Heston model, X will correspond 
to values of the stock price at discrete time intervals along a path, and Y will be a set of 
state variables recording information about the variance path. To derive the LR and the 
PW estimator, Broadie and Kaya differentiate inside the expectation operator in (0). If the 
interchange of differentiation and integration is justified for the left hand side of (0, it is 
also justified for the right hand side. Though the PW approach can be applied to the left 
hand side of (0, in ifTJl it is applied to the right hand side, to ensure the computational 
times are comparable. 

To show how to apply the approach, we fix a partition of the time interval, = to < 
t\ < ■ ■ ■ < t c j = T . As in the previous sections of the paper, we are interested in pricing 
a path-dependent option whose payoff is a function of the stock price vector (S t0 , . ..,S td ). 
Assuming that we have simulated a path of the variance process using the first two steps 
of the algorithms in Sections [5] and [6] and consider two consecutive times f; < ti. From 
Sections[5]and[6l we obtain values of f t j V$ds and f t j \fV~ s dW}. Define 



{l-p 2 )f t >V s ds 



and 



:= exp ( f ' V s ds + p f 1 VV s dW s l 



2 

Then given S tj , and the variance path, the value of St- can be expressed as 

2 



S t] = S t , $j exp ( (r - ^ ) (tj -tt) + Oj ^Tf^TiZ 



where Z is a standard normal random variable. Hence we take Y in © to be the variance 
path, and hence reduce the distribution in the inner expectation to a sequence of lognormal 
random variables, which is crucial to the LR method, as we need to be able to compute the 
score function. Again, as an example of a payoff, we consider again an Asian option with 
strike K, expiry T and payoff Lf=i ~ ^) + ■ Once the path wise and likelihood ratio 
estimates are given, the way to employ the algorithm from Section|5]becomes obvious. As 
such, to simplify notation, we set S = X Y?i=i St r define the time increment as Af,- := f, — f,_ i , 
and pose 



d: = 



(log(5,/(5,_ 1 ^,))-('-3^) Ar ' 



where af is the variance between f,_i and f,. Then the pathwise and likelihood ratio esti- 
mates are as follows. 

7.1.1. Pathwise ( PW) estimators. 

Delta : e^l^A 



1 " 



Rho : e- rT l^ K \ -ZSrfi-nS-K) 
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7.1.2. Likelihood ratio estimators. 



Gamma 



Delta 




e- rT (S 




Rho 




7.2. Barrier options. Consider a barrier option with monitoring dates (t\ , ti , . . . , t%), 
where we choose h := 2™ for simplicity. If the option is not knocked out at maturity 
T = th, it pays off /(St)- Following the approach of Glasserman and Staum [26 1, because 
of the knock-out feature, (S t ) t >o takes a value in R + U A, where A is an absorbing state. If 
S crosses the barrier at time f,, the option is knocked out, and for all j > i, Sj = A. De- 
fine A, to be the indicator function l(S tj ^ A), so A, = 1 means that the option is alive 
at time f,-. Assume the barrier is a price level H < So, so A; = 1 if the stock price has not 
crossed below the barrier H by step i. A down-and-out call in this model has the discounted 
pay-off A m e~ rT (St — K) + , where K is the strike price, r the constant interest rate, T = f/, 
denotes maturity. This allows us to propose a naive approach to pricing a barrier option in 
Algorithm[8] 

Sampling conditionally on one-step survival uses Algorithm [8] but in step 3 we use © 
where U = ( 1 — p (S t .) ) + Vp (S^.), with V uniformly distributed on (0,1), and 



We refer the reader to Section 2.2 in J26 | and in particular equations (1 1), (12), and (13). 



Algorithm 8 Bridge version of the Broadie-Kaya exact simulation algorithm 
1: Simulate (V;,)'^! as in Algorithm[T]or Algorithmic 



where m(u,t) := r(t-u) - \ f u V s ds + p f u «/V s dW} and C 2 (u,t) := (1 -p 2 ) J' u V s ds, 
where U is uniformly distributed on (0, 1) and <t> denotes the standard normal cdf. 



7.3. Multi-asset stochastic volatility models. In this subsection, we show that the ap- 
proach discussed in this paper can also be extended to a multi-asset stochastic volatility 
setting. When studying multidimensional stochastic volatility models, it is convenient to 
specify the model in such a way that the resulting multidimensional stochastic volatility 
process is an affine process [20, 22 1. For affine stochastic volatility models, characteristic 
functions are known to satisfy a particular set of Ricatti equations, resulting in a tractable 
model in which path-independent European options can be efficiently priced using Fourier 





3: Simulate (S,,.)f =1 as follows: Given S tj , V tj , V t}+1 , and V s ds, set 




(6) 
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transforms [ 16 1. As we now demonstrate, the methodology developed in this paper is par- 
ticularly amenable to the affine structure of the model: To be precise, we introduce two 
stock price processes, S l and S 2 , and three variance processes, V 1 , V 2 , and V 3 . We model 
the covariation between the stock prices via the variance processes as 

dS x t = S] [rdt + \Jvj-dZ} + \/vjdzf\ , 

dSj = S 2 [rdt + \JvjdZ 2 + y/vjdZfj , 
where for i = 1 , 2, 3 we have 

dv' = Ki (6, - v t ) dt + a.^Jv/dw;, z\ = piw; + ^\-p 2 b\, 

where (W 1 ,W 2 ,W 3 ,B l ,B 2 ,B 3 ) is a standard six-dimensional Brownian motion. We point 
out that each stock price process corresponds to a Bi-Heston model, as introduced by 
Christoffersen, Heston, and Jacobs in ifTTl . The stochastic covariation between S l and 
S 2 is given by (5 1 ,5 2 ) r = f 'SlS 2 V^ds. Using the terminology of GDI , this model is an 
A3(5). As in Section[2] we use the following representation for stock prices: 




where i = 1,2. From this representation, it is now clear how to produce an algorithm which 
allows us to handle multiasset stochastic volatility models. We propose this as Algorithmic 



Algorithm 9 Exact Simulation Algorithm for the Multiasset model 

1: Simulate V t \V 2 ,V t 3 using the noncentral ^-distribution 
2: Simulate j* V}ds\V} , & V 2 ds\V 2 , j* V ?ds\V t 3 
3: For i = 1,2,3; Compute 

f JvjdW* = — (V/ - V l + [ K.V'ds - KM) 

JO v O"; Jo 

4: Simulate log(S/) ~ N(Hi(t),a 2 (t)), where 

Mi (t) = log(5 ) + rt-\f V}ds -\f V^ds + pi /' JvjdW} + p 3 /' V^^ 3 
2 JO 2 Jo Jo v JO v 

af (0 - (i-pf) fv}ds+ (i -pi) fvfa 

Jo Jo 
and log {S 2 )~N{pL 2 {t) 1 d 2 2 {t)), where 

M2 (0 = log(Sg) + rt- l -[ V 2 ds -\f V^ds + p 2 f JvjdW 2 + P 3 f JvjdW* 
2 Jo 2 Jo JO v JO v 

6 2 2 {t) = (\-p 2 ) ['v 2 ds+(l-p 2 ) fvfds. 
Jo Jo 
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We point out that Algorithm [9] samples the joint distribution (Sj ,S 2 ). However, from 
Algorithm [5] it is clear how to modify Algorithm [9] to allow for path-dependent payoffs, 
i.e., how to sample the joint distribution 

(sl l ,...,sj d ,sf l ,...,sj d ). 

7.4. 3/2 Model. The 3/2 model was introduced by Heston in (29), see also ifBl 1331 l39l . 
Under the 3/2 model, the stock price is given by 

dS t 



rdt + VVp dWf+y/VtJl-p 2 dW t 2 , 

Of 

dV t = KV t (6-Vt)dt + eV t 3/2 dW t l . 

Recently, Baldeaux showed how to modify the Broadie and Kaya approach from Sec- 
tion [2] to handle the 3/2 model: First, we introduce the process X t =V t -\t>Q, whose 
dynamics are given by 

dX, = (k + e 2 - kOX,) dt - ey%dW t l . (7) 
It is now easily verified that the stock price S u is given by 



S,exp(r(K-f)-£jf {X s )- l ds + pJ t (y/X?) ' dW} + y/l - p 2 j 1 dW 2 

Equation (0 shows thatX = {X t , t > 0} is a square-root process. We recall the algorithm 
presented in J3) for the exact simulation of the 3/2 model in Algorithm [10] As for the 

Algorithm 10 Exact Simulation Algorithm for the 3/2 model 

1: Simulate X u \X t using the noncentral % 2 -distribution 
2: Simulate J t u jf\X t ,X u 

3: Recover J" (^%) ^ dW} from 



%) 1 dw} = l -(i og (^] + ( K + e l) r^- K e( U -t) 



e \ \X U J \ 2 J Jt X s 



4: Simulate S u given S t , J" (\/X^) dW} and j" (X s ) ds via 



log(S u )~N^log(S t ) + r(u-t)-^J\x s )- l ds + p J\VZ) l dW^a 2 {t,u) ) . 
where ^ 

a 2 (t,u) = (l -p 2 ) [ x; l ds. 



Heston model, Step 2 is the difficult step and, as in Section [2] one proceeds by computing 
the relevant Laplace transform [3|. We can use the bridge constructions for the square- 
root process and the Brownian bridge from Section [5] for steps 1 and 3 of Algorithm [TOl 
respectively, and handle Step 2 as in 0. Consequently, we can use the techniques from 
Section [5] to also handle the 3/2 model for path-dependent options. 
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